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The detrending moving average (DMA) algorithm is a widely used technique to quantify the 
long-term correlations of non-stationary time series and the long-range correlations of fractal sur- 
faces, which contains a parameter 9 determining the position of the detrending window. We develop 
multifractal detrending moving average (MFDMA) algorithms for the analysis of one-dimensional 
multifractal measures and higher-dimensional multifractals, which is a generalization of the DMA 
method. The performance of the one-dimensional and two-dimensional MFDMA methods is investi- 
gated using synthetic multifractal measures with analytical solutions for backward (8 = 0), centered 
(8 = 0.5), and forward (8 — 1) detrending windows. We find that the estimated multifractal scaling 
exponent r(q) and the singularity spectrum f(a) are in good agreement with the theoretical values. 
In addition, the backward MFDMA method has the best performance, which provides the most 
accurate estimates of the scaling exponents with lowest error bars, while the centered MFDMA 
method has the worse performance. It is found that the backward MFDMA algorithm also out- 
performs the multifractal detrended fluctuation analysis (MFDFA). The one-dimensional backward 
MFDMA method is applied to analyzing the time series of Shanghai Stock Exchange Composite 
Index and its multifractal nature is confirmed. 

PACS numbers: 05.45. Df, 05.40.-a, 05.10.-a, 89.75. Da 



I. INTRODUCTION 

Fractals and multifractals are ubiquitous in natural 
and social sciences [H-Q- There are a large number of 
methods developed to characterize the properties of frac- 
tals and multifractals. The classic method is the Hurst 
analysis or rcscalcd range analysis (R/S) for time se- 
ries [J, [f| and fractal surfaces [f|. The wavelet trans- 
form module maxima (WTMM) method is a more pow- 
erful tool to address the multifractality [fl-tllj. even for 
high-dimensional multifractal measures in the fields of 
image technology and three-dimensional turbulence [12l — 
fl6j . Another popular method is the detrended fluctua- 
tion analysis (DFA), which has the advantages of easy 
implementation and robust estimation even for short sig- 
nals [17H19J . The DFA method was originally invented to 
study the long-range dependence in coding and noncod- 
ing DNA nucleotides sequence 12011 and then applied to 
time series in various fields |2lT[24| . The DFA algorithm 
was extended to analyze the multifractal time series, 
which is termed as multifractal detrended fluctuation 
analysis (MFDFA) [Hj]. These DFA and MFDFA meth- 
ods were also generalized to analyze high-dimensional 
fractals and multifractals [2a ]. 

A more recent method is based on the moving aver- 
age (MA) or mobile average technique 27[, which was 
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first proposed by Vandcwalle and Ausloos to estimate 
the Hurst exponent of self- affinity signals [28[ and further 
developed to the detrending moving average (DMA) by 
considering the second-order difference between the orig- 
inal signal and its moving average function [23 ] . Because 
the DMA method can be easily implemented to estimate 
the correlation properties of non-stationary series with- 
out any assumption, it is widely applied to the analysis of 
real-world time series [30rl37| and synthetic signals (38l — 
|40| . Recently, Carbone extended the one-dimensional 
DMA method to higher dimensions to estimate the Hurst 
exponents of higher-dimensional fractals [4l|, [42( . Exten- 
sive numerical experiments unveil that the performance 
of the DMA method are comparable to the DFA method 
with slightly different priorities under different situations 

mm. 

In this paper, we extend the DMA method to mul- 
tifractal detrending moving average (MFDMA), which 
is designed to analyze multifractal time series and 
multifractal surfaces. Further extensions to higher- 
dimensional versions are straightforward. The perfor- 
mance of the MFDMA algorithms is investigated us- 
ing synthetic multifractal measures with known multi- 
fractal properties. We also compare the performance of 
MFDMA with MFDFA, and find that MFDMA is supe- 
rior to MFDFA for multifractal analysis. 

The paper is organized as follows. In Sec. [H] we de- 
scribe the algorithm of one-dimensional MFDMA and 
show the results of numerical simulations. We also ap- 
ply the one-dimensional MFDMA to analyze the time 



series of intraday Shanghai Stock Exchange Composite 
Index (SSEC). In Sec. IIII1 we describe the algorithm of 
two-dimensional MFDMA and report the results of nu- 
merical simulations as well. We discuss and conclude in 
Sec. EH 



II. ONE-DIMENSIONAL MULTIFRACTAL 
DETRENDING MOVING AVERAGE ANALYSIS 

A. Algorithm 



Step 1. Consider a time series x(t), t = 1,2, 
We construct the sequence of cumulative sums 



!/(*) = $>(*), * = 1,2, 
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,N. 



(1) 



; = 1 



Step 2. Calculate the moving average function y(t) in 
a moving window [35| . 

r(n-i)(i-e)i 

y(t) = - E ^- fc )' ( 2 ) 

n * — ' 
k=-L(n-i)ej 

where n is the window size, [x\ is the largest integer not 
greater than x, \x~\ is the smallest integer not smaller 
than x, and 9 is the position parameter with the value 
varying in the range [0,1]. Hence, the moving average 
function considers \(n — 1)(1 — 9)~\ data points in the 
past and [(n — 1)9 \ points in the future. We consider 
three special cases in this paper. The first case 9 = 
refers to the backward moving average [39( , in which the 
moving average function y(t) is calculated over all the 
past n — 1 data points of the signal. The second case 
9 = 0.5 corresponds to the centered moving average j39{ . 
where y(t) contains half past and half future information 
in each window. The third case 9 = 1 is called the for- 
ward moving average, where y(t) considers the trend of 
n — 1 data points in the future. 

Step 3. Detrend the signal series by removing the mov- 
ing average function y(i) from y(i), and obtain the resid- 
ual sequence e(i) through 



e 00 = y(i)-y(i), 



(3) 



where n-[(n- 1)9] ^ i sC N - [(n - 1)0J . 

Step 4- The residual series e(i) is divided into N n 
disjoint segments with the same size n, where N n = 
[N/n — lj . Each segment can be denoted by e v such 
that e v (i) = e(l + i) for 1 ^ i ^ n, where I = (v — l)n. 
The root- mean-square function F v (n) with the segment 
size n can be calculated by 



i=\ 



!(0- 



(4) 



Step 5. The gth order overall fluctuation function 
F q (n) is determined as follows, 



F q (n) 




(5) 



where q can take any real value except for q = 0. When 
q = 0, we have 



1 JV„ 
]n[F (n)] = w ^2HFv(n)] 



(6) 



according to L'Hospital's rule. 

Step 6. Varying the values of segment size n, we can 
determine the power-law relation between the function 
F q (n) and the size scale n, which reads 



F q (n) ~ n hiq \ 



(7) 



According to the standard multifractal formalism, the 
multifractal scaling exponent r(q) can be used to char- 
acterize the multifractal nature, which reads 



r(q) = qh(q) - D f , 



(8) 



where Df is the fractal dimension of the geometric sup- 
port of the multifractal measure [25|. For time series 
analysis, we have Df = 1. If the scaling exponent func- 
tion r(q) is a nonlinear function of q, the signal has 
multifractal nature. It is easy to obtain the singular- 
ity strength function a(q) and the multifractal spectrum 
f(a) via the Legendre transform [4J] 



a(q) = dr(q)/dq 
f(q) = qa- r(q) 



(9) 



B. Numerical experiments 

In the numerical experiments, we generate one- 
dimensional multifractal measure to investigate the per- 
formance of MFDMA, which is compared with MFDFA. 
We apply the p-model, a multiplicative cascading pro- 
cess, to synthesize the multifractal measure [45||. Start- 
ing from a measure fj, uniformly distributed on an in- 
terval [0,1]. In the first step, the measure is redis- 
tributed on the interval, /Ui,i = yq>\ to the first half and 
/ii.2 = HP2 = m(1 — Pi) to the second half. One partitions 
it into two sub-lines with the same length. In the (k + 1)- 
th step, the measure Hk,i on each of the 2 fe line segments 
is redistributed into two parts, where fJ,k+i,2i-i = l^k,tPi 
and jik+i,2i = t i k,iP2- We repeat the procedure for 14 
times and finally generate the one-dimensional multi- 
fractal measure with the length 2 14 . In this paper, we 
present the results when the parameters are p\ = 0.3 
and P2 = 0.7. The results for other parameters are qual- 
itatively the same. 




FIG. 1. (Color online) Multifractal analysis of the one-dimensional multifractal binomial measure using the three MFDMA 
algorithms and the MFDFA approach, (a) Power-law dependence of the fluctuation functions F q (n) with respect to the scale 
n for q = —4, q — 0, and q = 4. The straight lines are the best power-law fits to the data. The results have been translated 
vertically for better visibility, (b) Multifractal mass exponents r(q) obtained from the MFDMA and MFDFA methods with 
the theoretical curve shown as a solid line, (c) Differences Ar(q) between the estimated mass exponents and their theoretical 
values for the four algorithms, (d) Multifractal spectra f(a) with respect to the singularity strength a for the four methods. 
The continuous curve is the theoretical multifractal spectrum. 



We calculate the fluctuation function F q (n) of the syn- 
thetical multifractal measure using the MFDMA and 
MFDFA methods, and present the fluctuation function 
F q (n) in Fig. Q^a). We find that the function F q {n) well 
scales with the scale size n. Using the least squares fitting 
method, we obtain the slopes h(q) for MFDMA (9 = 0, 
= 0.5 and = 1) and MFDFA respectively, which are 
illustrated in Table |U It is found that the error bars of 
the three MFDMA algorithms are all smaller than the 
MFDFA method, which implies that it is easier to deter- 
mine the scaling ranges for the MFDMA algorithms. In 
most cases, the algorithms underestimate the h(q) val- 
ues and the backward MFDMA approach gives the best 
estimates. There is an interesting feature in Fig. [TJa) 
showing evident log-periodic oscillations in the MFDFA 
Fq{n) curves, which is intrinsic for the multifractal bino- 
mial measure |46( . 

We plot the multifractal scaling exponents r(q) ob- 
tained from MFDMA (6 = 0, 6 = 0.5 and = 1) and 



TABLE I. The MFDMA exponents h(q) for q = -4, -2, 0, 
2, and 4 of the one-dimensional synthetic multifractal mea- 
sure with the parameters p\ = 0.3 and P2 — 0.7 using the 
MFDMA (0 = 0,0 = 0.5 and = 1) and MFDFA methods. 
The numbers in the parentheses are the standard errors of 
the regression coefficient estimates using the t-test at the 5% 
significance level. 



MFDMA 



= 



= 0.5 



= 1 



MFDFA Analytic 



-4 1.505(4) 1.401(12) 1.496(2) 

-2 1.354(3) 1.249( 8) 1.337(4) 

1.114(4) 1.022( 5) 1.096(5) 

2 0.874(6) 0.788( 3) 0.859(5) 

4 0.749(9) 0.667( 4) 0.736(6) 



1.490(17) 1.499 

1.326( 9) 1.359 

1.074( 6) 1.126 

0.804(11) 0.893 

0.670(15) 0.753 



MFDFA in Fig.^b). The theoretical formula of r(q) of 
the multifractal measure generated by the p-model dis- 
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FIG. 2. (Color online) Multifractal analysis of the 5-min return time series of the SSEC index using the backward MFDMA 
method, (a) Power-law dependence of the fluctuation functions F(n) with respect to the scale n. The solid lines are the 
least-squares fits to the data. The results corresponding to q — —2, q — 0, q — 2 and q = 4 have been translated vertically 
for clarity, (b) Multifractal spectra f(a) of the raw return series of SSEC and its shuffled series. Inset: Multifractal scaling 
exponents r(q) as a function of q. 



cussed above can be expressed by [44 



7th (g) 



HpI + pD 

In 2 



(10) 



which has been illustrated in Fig.QJb) as well. In order to 
quantitatively evaluate the performance of MFDMA and 
MFDFA, we calculate the relative estimation errors of the 
numerical values of r(q) in reference to the corresponding 
theoretical values 7th (?) 



Ar(q) =r(q) -T th (q), 



(11) 



which arc shown in Fig. QJc). When < q ^ 4, all the 
four methods underestimate the r(q) exponents. In con- 
trast, when —4 ^ q < 0, these methods overestimate the 
scaling exponents with a few exceptions for the backward 
MFDMA case. It is evident that the backward MFDMA 
method (0 = 0) gives the most accurate estimation of 
the exponents, the forward MFDMA method (9 = 1) has 
the second best performance, and the centered MFDMA 
method (9 = 0.5) performs worst. We stress that both 
the backward and the forward MFDMA methods outper- 
form the MFDFA approach. 

According to the Legendre transform, we can numeri- 
cally calculate the singularity strength functions a(q) and 
the multifractal spectrum functions f(a) with Eq. @ 
for the four methods, which are depicted in Fig. [TJd). 
We also show the theoretical singularity spectrum as a 
continuous curve for comparison, where the singularity 
strength function a(q) can be calculated as follows 



a t h{q) 



p\\npi +pl\np 2 
{p\+p q 2 )\u2 



and the multifractal spectrum is 



/th(<?) = - 



qp\ lnpi + qp\ lnp 2 . ln(p? + p q 2 ) 



{p\+pl)\n2 



In 2 



(12) 



(13) 



It also shows that the order of algorithm performance 
is the following: backward MFDMA, forward MFDMA, 
MFDFA, and centered MFDMA. 



C. Application to intraday SSEC time series 

We now apply the one-dimensional backward MFDMA 
method to the study of the multifractal properties of 5- 
min return series of Shanghai Stock Exchange Composite 
Index (SSEC) within the time period from January 2003 
to April 2008. The number of data points is about 10 5 . In 
Fig-Uta), we present the fluctuation function F q (n) with 
respect to the scale n. It is found that the function F q {n) 
and the scale n have a sound power-law relation. The 
MFDMA exponents h{q) can be estimated by the slopes 
of the straight lines illustrated in Fig.^a) for different q. 
Using the least-squares fitting method, we have h(— 4) = 
0.660±0.005, h(-2) = 0.624±0.002, h(0) = 0.591±0.001, 
h{2) = 0.531±0.003, and h(4) = 0.493±0.008. Since h(2) 
is close to 0.5, it is confirmed that the return time series 
of the SSEC is almost uncorrected, which is consistent 
with the earlier empirical results. 

Figure[2]Jb) shows the multifractal spectrum /(a) of re- 
turn series of SSEC. The strength of mult ifr act ality can 
be characterized by the span of the multifractal singular- 
ity strength function, that is, Aa = cv max — a m i„. If Act is 
large, the return series owns multifractality, while the re- 
turn series is almost monofractal if Aa gets close to zero. 
We observe in the figure that Aa = 0.72 - 0.39 = 0.33, 
which means that the return series of SSEC has mul- 
tifractal nature. The nonlinearity of the r(q) curve in 
the inset confirms the presence of multifractality in the 
return series. 

We shuffle the return time series and reconstruct a 



surrogate index. We then perform backward MFDMA 
analysis of the surrogate index. The results are also il- 
lustrated in Fig. EJJb). We find that the singularity spec- 
trum shrinks much and the r(q) function is almost linear, 
which implies that the fat-tailedness of the returns plays 
a minor role in producing the observed multifractality 
and the correlation structure is the main cause of mul- 
tifractality. This observation is quite interesting since it 
is different from the conclusion that the fat-tailedness of 
the returns is the main origin of multifractality according 
to the MFDFA method Eg. 



III. TWO-DIMENSIONAL MULTIFRACTAL 
DETRENDING MOVING AVERAGE ANALYSIS 

A. Algorithm 

The two-dimensional MFDMA algorithm is used to 
investigate possible multifractal properties of surfaces, 
which can be denoted by a two-dimensional matrix 
X(h,i 2 ) with h = 1,2,--- ,Ni and i 2 = 1,2,- •• ,N 2 . 
The algorithm is described as follows. 

Step 1. Calculate the sum Y(ii, i 2 ) in a sliding window 
with size n\ x n 2 , where n\ ^ i\ ^ N\ — [( n i — l)0ij 
and n 2 ^ i 2 ^ N 2 — [{n 2 — 1)0 2 J. The two position 
parameters 9\ and 9 2 vary in the range [0,1]. Specifi- 
cally, we extract a sub-matrix Z(u±, u 2 ) with size n\ x n 2 
from the matrix X, where i% — n% + 1 ^ u\ ^ i\ and 
i 2 — n 2 + 1 ^ u 2 ^ i 2 . We can calculate the sum Y(i\, i 2 ) 
of Z as follows, 






(14) 



Step 2. Determine the moving average function 
Y(i\,i 2 ), where m ^ i\ ^ N\ — \_(ni — l)9i\ and n 2 ^ 
i 2 ^ N 2 — [(n 2 — 1)6*2] ■ We first extract a sub-matrix 
W(ki,k 2 ) with size ri\ x n 2 from the matrix X, where 
h - \(ni - 1)(1 - 6>i)l sC fci < fci + [(ni - l)0i J and 
k 2 - \(n 2 - 1)(1 - 2 )] < fc 2 < fc 2 _f L(n 2 - 1)0 2 J. Then 

we calculate the cumulative sum W(mi,m 2 ) of W, 



W(mi, m 2 ) 



mi ran 

EE 

di=ld 2 =l 



W(di,d2), 



(15) 



where 1 ^ mi Sj ni and 1 $J m 2 ^ n 2 . The moving 
average function Y(ii, i 2 ) can be calculated as follows, 



Y(i u i 2 ) = 



1 



nin 2 



ni «2 

E 

mi-l 7712 — 1 



^ W(mi,m 2 ). (16) 



££ep 5. Detrend the matrix by removing the moving 
average function Y(ii,i 2 ) from Y(i±,i 2 ), and obtain the 
residual matrix e(ii,i 2 ) as follows, 



where n\ ^ i\ ^ iVi — [(^l — l)0ij and n 2 ^ i 2 $J 

A r 2 " L(«2 " 1)02 J ; 

S'iep 4- The residual matrix e(ii,i 2 ) is partitioned into 
N ni x N n2 disjoint rectangle segments of the same size 
n\ x n 2 , where N ni = [(Ni — ni(l + 0i))/nij and N ri2 = 
l(N 2 — n 2 (l + 9 2 ))/n 2 \. Each segment can be denoted 
by e Vl , V2 such that e Vl . v . 2 {ii,i 2 ) = e(h +ii,l 2 + i 2 ) for 
I ^ ii ^ ni and I ^ i 2 ^ n 2 , where l\ = (v\ — l)ni and 
h = {v 2 — l)n 2 . The detrended fluctuation F Vl , V2 {n\, n 2 ) 
of segment e UI)W2 (ii,i2) can be calculated as follows, 



F^>i,n 2 ) = — - ^ ^ < )VB (ii,i 2 ). (18) 

'^1 '''2 



-1 "T "2 

— EE 

,Tl n ^— ' -^— ' 
21 — 1 22 — 1 

Step 5. The gth order overall fluctuation function 
F q {n) is calculated as follows, 



F q (n) = 



N ni N n2 

EE 



N ni N n2 ^ ^ 

1 z v 1— 1 V2 — 1 



F^^^i' 7 ^) 



(19) 



where n 2 = 2( n i + n l) an d 9 can take any real values 
except for q = 0. When q = 0, we have 



ln[F (n)] = 



1 



N ni N n2 

X) EH^i.^Cni.na)], (20) 

1 772 ij 1 — X D2 — 1 



according to L'Hospital's rule. 

Step 6. Varying the segment sizes n\ and n 2 , we are 
able to determine the power-law relation between the 
fluctuation function F q {n) and the scale n, 



F q {n) 



Mi) 



(21) 



e(h,i 2 ) = Y{ix,i 2 ) - Y{ix,i 2 ), 



(17) 



In this paper, we particularly adopt n = ri\ = n 2 and 
9 = 9\ = 9 2 for the isotropic implementation of the tow- 
dimensional MFDMA algorithm. Applying Eqs. ([5]) and 
([9]). we can obtain the multifractal scaling exponent r(g), 
the singularity strength function a(q) and the multifrac- 
tal spectrum f(a), respectively. For the two-dimensional 
multifractal measures, we have Df = 2 in Eq. ([S]). 



B. Numerical experiments 

In order to investigate the performance of the two- 
dimensional MFDMA methods, we adopt the multiplica- 
tive cascading process to synthesize the two-dimensional 
multifractal measure. The process begins with a square, 
and we partition it into four sub-squares with the same 
size. We then assign four proportions of measure pi, p 2 , 
P3 and pi to them (s.t. p\ +p 2 +P3 -\-pi = 1). Each sub- 
square is further partitioned into four smaller squares 
and the measure is re-assigned with the same propor- 
tions. The procedure is repeated 10 times and we finally 
generate the two-dimensional multifractal measure with 
size 1024 x 1024. In the paper, the model parameters are 
p x = 0.1, p 2 = 0.2, p 3 = 0.3, and p 4 = 0.4. 
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FIG. 3. (Color online) Multifractal analysis of the two-dimensional multifractal measure using the three MFDMA algorithms 
and the MFDFA approach, (a) Power-law dependence of the fluctuation functions F q (n) with respect to the scale n for q = —4, 
q = 0, and q = 4. The straight lines are the best power-law fits to the data. The results have been translated vertically for 
better visibility, (b) Multifractal mass exponents r(q) obtained from the MFDMA and MFDFA methods with the theoretical 
curve shown as a solid line, (c) Differences Ar(q) between the estimated mass exponents and their theoretical values for the 
four algorithms, (d) Multifractal spectra f(a) with respect to the singularity strength a for the four methods. The continuous 
curve is the theoretical multifractal spectrum. 



In Fig. [3ja) , we depict the fluctuation functions F q (n) 
of the two-dimensional multifractal measure for three 
MFDMA algorit hms w ith 6 = 0, 6 = 0.5 and = 1 
described in Sec. IIII Al and the two-dimensional MFDFA 
method [26|]. We find that every F q (n) function scales 
excellently as a power law of n. Linear least-squares re- 
gressions of ln[F n (g)] against ln(n) for each q give the 
estimates of h(q). Table [TTI shows the resultant values of 
h(q) for q = —4, -2, 0, 2, and 4. Similar to the results 
of the one-dimensional MFDMA and MFDFA methods, 
we find that the three two-dimensional MFDMA algo- 
rithms give better power-law scaling relations than the 
MFDFA algorithm with smaller standard errors in the 
parentheses except for the forward MFDMA with 6 = 1 
when q = —4. Compared with the theoretical values in 
the last column of Table HH the backward MFDMA with 
6 = gives the most accurate estimates and thus has the 
best performance. 

In Fig. [3](b) , we plot the multifractal scaling exponents 



TABLE II. The MFDMA exponents h(q) for q = -4, -2, 0, 2, 
and 4 of the two-dimensional synthetic multifractal measure 
with the parameters pi = 0.1, pi = 0.2, pz = 0.3, andp4 = 0.4 
using the MFDMA (6 = 0, 6 = 0.5 and 6=1) and MFDFA 
methods. The numbers in the parentheses are the standard 
errors of the regression coefficient estimates using the t-test 
at the 5% significance level. 



MFDMA 







0.5 6 



MFDFA Analytic 



-4 2.850( 7) 2.857(11) 2.860(42) 

-2 2.534( 6) 2.505(14) 2.503(22) 

2.114(10) 2.041(13) 2.017(15) 

2 1.829(19) 1.751(10) 1.720(10) 

4 1.688(25) 1.615(10) 1.586(11) 



2.780(20) 


2.849 


2.461(23) 


2.577 


2.064(20) 


2.176 


1.769(30) 


1.869 


1.616(43) 


1.705 



r(q) estimated from the three MFDMA algorithms and 
the MFDFA method. The theoretical results are also 



shown for comparison, which are calculated as follows 



Tth(g) 



In 2 



(22) 



According to Fig- EJ^b), we find that the four r{q) curves 
estimated from the MFDMA and MFDFA methods arc 
also close to the theoretical values. 

To have a better visibility, we plot the difference func- 
tions Ar(q) = r(q) — r t h{q) in Fig. [3l[b). For positive 
values of q, the algorithms underestimate the r(q) val- 
ues. It is clear that the backward MFDMA with 8 = 
performs best, the forward MFDMA with 8 = 1 per- 
forms worst, and the MFDFA outperforms the centered 
MFDMA with 8 = 0.5. For negative values of q, most 
At values are larger than 0. We find that the backward 
MFDMA has the best performance, the MFDFA has the 
worse performance, and the performance of the centered 
and forward MFDMA methods are comparable to each 
other. These findings are further confirmed by the results 
illustrated in Fig. EJd) for the multifractal spectra. 



IV. DISCUSSION AND CONCLUSION 

In this paper, we have developed the detrending mov- 
ing average techniques to make them suitable for the 
analysis of multifractal measures in one and two dimen- 
sions. Extensions to higher dimensions are straightfor- 
ward. The performance of these MFDMA algorithms is 
tested based on numerical experiments of synthetic mul- 
tifractal measures with known theoretical multifractal 
properties generated according to multiplicative cascad- 
ing processes. We also present the results of MFDFA for 
comparison. Our main conclusion is that the backward 
MFDMA with the parameter 8 = exhibits the best per- 
formance when compared with the centered MFDMA, 
forward MFDMA, and MFDFA, because it gives better 
power-law scaling in the fluctuation functions and more 
accurate estimates of the multifractal scaling exponents 
and the singularity spectrum. 

For the one-dimensional MFDMA version, we find that 
MFDMA gives a more reliable regression in the log-log 
plot of the fluctuation function F q (n) with respect to 
the scale n than MFDFA. Comparing with the theoreti- 
cal formulas of r(q) and f(a), the backward MFDMA 



performs best, the centered MFDMA performs worst, 
and the forward MFDMA outperforms the MFDFA. The 
backward MFDMA with 8 = is applied to the analysis 
of the return time series of SSEC and confirms that the 
return series exhibits multifractal nature, which is not 
caused by the fat-tailedness of the return distribution. 

For the two-dimensional MFDMA version, we also 
find that MFDMA gives a more reliable regression than 
MFDFA when testing on the two-dimensional synthetic 
multifractal measure. The estimates of r(q) and f(a) 
well agree with the theoretical values for the backward 
MFDMA with = 0, which is the best estimator. The 
centered and forward MFDMA methods with 8 = 0.5 
and 8 = 1 give a better estimation than MFDFA in the 
negative range of q, while they are worse than MFDFA 
for positive q. 

Technically, it is crucial to emphasize that, the win- 
dow size (n for one-dimensional MFDMA or ri\ x n^ for 
two-dimensional MFDMA) used to determine the mov- 
ing average function (y(i) for one-dimensional version or 
Y(ii,i2) for two-dimensional version) in Step 2 must be 
identical to the partitioning segment size in Step 4. If the 
window size is not equal to the segment size, the fluctua- 
tion function F q (n) does not show power-law dependence 
on the scale n, and the multifractal scaling exponent r(q) 
and the multifractal spectrum f(a) estimated from the 
MFDMA remarkably deviate from the theoretical ones. 

Finally, we would like to stress that there are tremen- 
dous potential applications of the backward MFDMA 
method in multifractal analysis, due to the better per- 
formance of this algorithm compared with the exten- 
sively used MFDFA approach. Possible applications in- 
clude time series (one-dimensional), fracture surfaces, 
landscapes, clouds, and many other images possess- 
ing self-similar properties (two-dimensional), tempera- 
ture fields and concentration fields (three-dimensional), 
and strange attractors in nonlinear dynamics (higher- 
di mensional) . The M ATLA B codes are publicly available 
at http://arxiv.org/abs/1005.0877 for download. 
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Appendix: MATLAB codes for Multifractal Detrending Moving Average (MFDMA) analysis 

These MATLAB codes are not submitted for publication. They are available as the appendix of the manuscript 
on the arXiv. Rather than copying directly from the generated PDF file, the readers are advised to download the 
source files and copy the MATLAB codes from the tex file. Questions and suggestions, if any, should be addressed to 
gfgu@ecust.edu.cn or wxzhou@ecust.edu.cn. 

% GFGUMFDMAJD.m 



function [n ,Fq, tau , alpha , f] = GFGU_MFDMA_1D ( x , n.min , n_max ,N, thcta , q) 

% 

% The 'procedure [n , Fq , tau , alpha , f ]~GFGU-MFDMA_lD(x , n_min , n_max ,N, theta , q) is 

% used to calculate the multifractal properties of one— dimensional time series 

% 

% Input : 

% x: the time series we considered 

% n_min : the lower bound of the segment size n 

% njmax: the upper bound of the segment size n 

% N: the length of n, that is , the data points in the plot of Fq VS n 



% theta : the position parameter of the moving window 

% q: multijractal order 

% 

% Output: 

% n: segment size series 

% Fq: q—th order fluctuation function 

% tau: multifractal scaling exponent 

% alpha: multifractal singularity strength function 

% f: multifractal spectrum 

% 

% The procedure works as follows: 

% 1) Construct the cumulative sum y. 

% 2) For each n, calculate the moving average function \widetilde {y} . 

% 3) Determine the residual e by detrending \widetilde {y} from y. 

% 4) Estimate the root— mean— square function F. 

% 5) Calculate the q—th order overall fluctuation function Fq. 

% 6) Calculate the multifractal scaling exponent tau(q). 

% 7) Calculate the singularity strength function alpha(q) and spectrum f(alpha). 

% 

% Note: 

% 1) The window size and the segment size must be identical. 

% 2) The lower bound n_min would better be selected around 10. 

% 3) The upper bound n_max would better be s elect ed around 10% of the length of 

% time series . 

% 4) N would better be s eleceted in the range [20 ,40}. 

% 5) The parameter theta varies in the range [0,1]. Theta ~ corresponds to 

% backward MFDMA, and theta — 0.5 corresponds to the centered MFDMA, and 

% theta — 1 corresponds to the forward MFDMA. We recommend theta =0. 

% 

% Example : 

% [n,Fq,tau , alpha . f ]=GFGUMFDMAJD(x , 1 , round ( length (x) / 10) ,30,0, -4:0.1:4) ; 

% 

if size (x ,2 ) = 1 

x = x ' ; 
end 

M = length (x) ; 

MIN = loglO(n.min) ; 

MAX = loglO(n_max) ; 

n = ( unique (round (logspace (MM, MAX,N) ))) '; 

% Construct the cumulative sum y 
y = cumsum(x) ; 

for i = 1: length (n) 

lgth = n(i ,1); 

% Calculate the moving average function \widetilde {y} 
yl = zeros (1 ,M-lgth + l); 
for j = l:M-lgth+l 

yl(j) = mean(y(j : j+lgth-1)) ; 
end 

% Determine the residual e 

c=y(max(l , floor (lgth *(1- theta ) ) ) :max(l , floor ( lgth *(1- theta )) )+length (yl ) -l)-yl 

% Estimate the root— mean— square function F 
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for k = l: floor ( length ( c )/ lgth ) 

F{i}(k)=sqrt(mean(e((k-l)*lgth + l:k*lgth) ."2)) 
end 
end 



% Calculate the q—th order overall fluctuation function Fq 
for i = l:length (q) 

for j = l:length (F) 

f=F{J}; 

if q(i) = 

Fq( j , i )=exp (0. 5* mean ( log( f . "2) )) ; 
else 

Fq(j ,i)=(mean(f."q(i)))"(l/q(i)); 
end 
end 
end 



% Calculate the multifractal scaling exponent tau(q) 
for i = 1 : s i z e ( Fq , 2 ) 

fq=Fq(: ,i); 

r=rcgstats(log(fq) , log(n) , 'linear ' ,{ ' tstat '}) ; 

k=r . tstat . beta ( 2 ) ; 

h(i,l)=k; 
end 
tau=h. *q' — 1; 

% Calculate the singularity strength function alpha(q) and spectrum f(alpha) 

dx = 7; 

dx=fix((dx-l)/2) ; 

for i=dx + l:length ( tau )— dx 

xx=q ( i — dx : i+dx ) ; 

yy=tau ( i— dx : i+dx) ; 

r=regst ats (yy,xx, 'linear ' ,{ 'tstat '}) ; 

alpha(i , l)=r . tstat . beta (2) ; 
end 

alpha=alpha (dx+l:end) ; 
f=q(dx + l:end— dx) ' . * alpha— tau (dx + l:end— dx) ; 

% GFGUMFDMAJiD .m 

function [n ,Fq, tau , alpha , f] = GFGUJMFDMA_2D(X, n_min ,n_max ,N, theta , q) 

% 

% The procedure [n , Fq , tau , alpha , f ' ] —GFGU-MFDMAJiD ( x , n_min , njmax ,N, theta , q) 

% is used to calculate the multifractal properties of two-dimensional 

% multifractal measures . 

% 

% Inp u t : 

% X: the two— dimensional multifractal measures we considered . 

% n_min : the lower bound of the segment size n 

% njmax: the upper bound of the segment size n 

% N: the length of n, that is , the data points in the plot of Fq VS n 

% theta: the position parameter of the moving window 

% q: multifractal order 

% 
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% Output: 

% n: segment size series 

% Fq: q—th order fluctuation function 

% tau: multifractal scaling exponent 

% alpha: multifractal singularity strength function 

% f: multifractal spectrum 

% 

% The procedure works as follows: 

% 1) For each n, construct the cumulative sum Y in a moving window. 

% 2) Calculate the moving average function \widetilde {Y} . 

% 3) Determine the residual e by detrending \widetilde {Y\ from Y. 

% 4) Estimate the root— mean— square function F. 

% 5) Calculate the q—th order overall fluctuation function Fq. 

% 6) Calculate the multifractal scaling exponent tau(q). 

% 7) Calculate the singularity strength function alpha(q) and spectrum f(alpha). 

% 

% Note: 

% 1) The window size and the segment size must be identical. 

% 2) The lower bound n_min would better be selected around 10. 

% 3) The upper bound n_max would better be selected around 10% of min( size (X) ) . 

% 4) N would better be seleceted in the range [20 ,40]. 

% 5) The parameter theta varies in the range [0,1] , and we have 

% theta=theta_l—theta_2. Theta — corresponds to backward MFDMA, and 

% theta = 0.5 corresponds to the centered MFDMA, and theta = 1 corresponds to 

% the forward MFDMA. We recommend theta —0. 

% 6) In the procedure , we have n=n_l=n_2 for the segment size. 

% 

% Example : 

% [n,Fq,tau , alpha , f ]=GFGUJMFDMA^D(X, 1 , round (min( size (X) )/ 10) ,30,0, -4:0.1:4) ; 

% 

Nl=size(X,l) ; 

N2=size(X,2) ; 

MM=loglO(n_min) ; 

MAX=loglO(n_max) ; 

n=(unique (round (logspace (MM, MAX,N) )) ) '; 

for i =l:length (n) 
lgth=n(i ,1) ; 

Y= zeros (Nl-lgth+l,N2-lgth+l) ; 
Yl = zeros (Nl-lgth + l,N2-lgth + l) ; 
for j = l:Nl-lgth+l 

for k = l:N2-lgth+l 

Z =X(j:j+lgth-l, k:k+lgth-l); 

Zl = (cumsum( (cumsum(Z) ) ')) '; 

% Construct the cumulative sum Y 
Y(j ,k)=Zl(end,end) ; 

% Calculate the moving average function \widetilde {Y} 
Yl(j ,k)=mean(Zl(:) ) ; 
end 
end 

% Determine the residual e 
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xO = l:size(Y,l )— min( floor(lgth*theta) , lgth — 1) ; 
yO = l: size (Y, 2 )-min( floor (lgth*theta) , lgth-1); 
xl=s iz e ( Yl,l) -length (xO) +1: size (Yl , 1 ) ; 
yl=size(Yl,2)-length(yO)+l:size(Yl,2) ; 
e=Y(xO,yO)-Yl( X l,yl); 

% Estimate the root— mean—square function F 
for kl = l:floor(size(e,l)/lgth) 

for k2 = l:floor(size(e,2)/lgth) 

E=e((kl-l)*lgth + l:kl*lgth , (k2-l)*lgth+l:k2*lgth ) ; 
F{i}(kl ,k2)=sqrt(mean(E(:) ."2) ) ; 
end 
end 
end 



% Calculate the q—th order overall fluctuation function Fq 
for i =l:length (q) 

for j = l:length(F) 

f=F{j }(:) ; 
if q(i) = 

Fq( j , i )=exp (0. 5* mean ( log(f."2))); 
else 

Fq(j ,i)=(mean(f."q(i)))"(l/q(i)); 
end 
end 
end 



% Calculate the multifractal scaling exponent tau(q) 
for i =l:size (Fq,2 ) 

fq=Fq(:,i); 

r=regstats (log(fq) , log(n) , 'linear ' ,{ ' tstat '}) ; 

k=r . tstat . beta (2) ; 

h(i,l)=k; 
end 
tau=h. *q' —2; 

% Calculate the singularity strength function alpha(q) and spectrum f(alpha) 

dx = 7; 

dx=fix((dx-l)/2) ; 

for i=dx + l:length ( tau )— dx 

xx=q ( i — dx : i+dx ) ; 

yy=t au ( i — dx : i +dx ) ; 

r=regstats(yy,xx. 'linear ' ,{ 'tstat '}) ; 

alpha(i ,l)=r. tstat . beta (2) ; 
end 

alpha=alpha (dx+l:end) ; 
f=q(dx + l:end— dx) '.* alpha— tau (dx + l:end— dx) ; 

I 



